# coding=gbk
"""
Created on 2016/6/23 0023

@author: Blunce

@about:

"""
from itertools import count, islice
import matplotlib.pyplot as plt
import numpy as np
import pandas as pd
from numpy.polynomial.polynomial import polyval, polyfit
from utils.BinomialMultiplicativeCascades import BMC


class MFDFA:
    def __init__(self, x):
        self.__x = x
        self.__m = 2
        self.__scales = np.arange(10, len(x) / 5, (len(x) / 5 - 10) / 100)
        self.__qs = np.arange(-10, 10.5, 0.5)
        self.__fq = np.zeros((len(self.__scales), len(self.__qs)))
        self.__hq = np.zeros(len(self.__qs))
        self.__tq = np.zeros(len(self.__qs))
        self.__a = np.zeros(len(self.__qs) - 1)
        self.__f = np.zeros(len(self.__qs) - 1)

    def setScales(self, scales):
        self.__scales = scales
        self.__fq = np.zeros((len(self.__scales), len(self.__qs)))

    def setM(self, m):
        self.__m = m

    def setQs(self, qs):
        self.__qs = qs
        self.__fq = np.zeros((len(self.__scales), len(self.__qs)))
        self.__hq = np.zeros(len(self.__qs))
        self.__tq = np.zeros(len(self.__qs))
        self.__a = np.zeros(len(self.__qs) - 1)
        self.__f = np.zeros(len(self.__qs) - 1)

    def process(self):
        x_cum = np.cumsum(self.__x - np.mean(self.__x))
        self.__rms(x_cum)
        lnFq = np.log(self.__fq)
        lnScales = np.log(self.__scales)
        for qi in range(len(self.__qs)):
            line_fit = polyfit(lnScales, lnFq[:, qi], 1)
            self.__hq[qi] = line_fit[1]
        self.__tq = self.__hq * self.__qs - 1
        self.__a = np.diff(self.__tq) / (self.__qs[2] - self.__qs[1])
        self.__f = self.__a * self.__qs[:-1] - self.__tq[:-1]
        return self.__hq, self.__tq, self.__a, self.__f

    def __rms(self, cumsum):
        cumsum_index = np.arange(len(cumsum))
        for si, scale in enumerate(self.__scales):
            ns = int(np.ceil(len(self.__x) / float(scale)))
            rest = len(self.__x) % scale
            div = len(self.__x) / scale
            l_rms = np.zeros(ns * 2)
            for i in range(ns * 2):
                tmp, tmp_index = [], []
                if i == ns - 1:
                    if rest != 0:
                        tmp = cumsum[div * scale:]
                        tmp_index = cumsum_index[div * scale:]
                    else:
                        tmp = cumsum[(div - 1) * scale:]
                        tmp_index = cumsum_index[(div - 1) * scale:]
                elif i == ns:
                    if rest != 0:
                        tmp = cumsum[:rest]
                        tmp_index = cumsum_index[:rest]
                    else:
                        tmp = cumsum[:scale]
                        tmp_index = cumsum_index[:scale]
                elif ns - 1 < i <= 2 * ns - 1:
                    if rest != 0:
                        tmp = cumsum[
                              rest + scale * (i - ns - 1):rest + scale * (
                                  i - ns)]
                        tmp_index = cumsum_index[
                                    rest + scale * (i - ns - 1):rest + scale * (
                                        i - ns)]
                    else:
                        tmp = cumsum[
                              rest + scale * (i - ns):rest + scale * (
                                  i - ns + 1)]
                        tmp_index = cumsum_index[
                                    rest + scale * (i - ns):rest + scale * (
                                        i - ns + 1)]
                elif 0 <= i < ns - 1:
                    tmp = cumsum[i * scale:(i + 1) * scale]
                    tmp_index = cumsum_index[i * scale:(i + 1) * scale]
                line_fit = polyfit(tmp_index, tmp, self.__m)
                tmp_fit = polyval(tmp_index, line_fit)
                l_rms[i] = np.mean((tmp - tmp_fit) ** 2)
            for qi, q in enumerate(self.__qs):
                if q == 0:
                    self.__fq[si, qi] = np.exp(np.mean(np.log(l_rms)) / 2.0)
                else:
                    self.__fq[si, qi] = (np.mean(l_rms ** (q / 2.0))) \
                                        ** (1.0 / q)
        return self.__fq

    def drawHq(self):
        plt.title("General Hurst Exponent")
        plt.plot(self.__qs, self.__hq)
        plt.show()

    def drawTq(self):
        plt.title("Mass Exponent")
        plt.plot(self.__qs, self.__tq)
        plt.show()

    def drawSpectrum(self):
        plt.title("Spectrum")
        plt.plot(self.__a, self.__f)
        plt.show()


def verify():
    bmc = BMC()
    x = bmc.process()
    a_theoretical, f_theoretical = bmc.getSpectrum()
    mfdfa = MFDFA(x)
    mfdfa.setM(1)
    mfdfa.setScales(np.arange(10, 5000, 50))
    hq_mfdfa, tq_mfdfa, a_mfdfa, f_mfdfa = mfdfa.process()
    plt.title("Verification of MFDFA")
    plt.plot(a_theoretical, f_theoretical, label="theoretical")
    plt.plot(a_mfdfa, f_mfdfa, label="MFDFA")
    plt.ylim(np.min([f_theoretical.min(), f_mfdfa.min()]) - 0.2,
             np.max([f_theoretical.max(), f_mfdfa.max()]) + 0.2)
    plt.xlim(np.min([a_theoretical.min(), a_mfdfa.min()]) - 0.2,
             np.max([a_theoretical.min(), a_mfdfa.max()]) + 0.2)
    plt.legend(loc="best")
    plt.show()


def demo():
    t = np.arange(0, 7, 0.01)
    signal = np.sin(t)
    mfdfa = MFDFA(signal)
    mfdfa.process()
    mfdfa.drawHq()
    mfdfa.drawTq()
    mfdfa.drawSpectrum()


def expmt_div(n=8, s=3):
    x = np.arange(n)
    print x
    scale = s
    ns = int(np.ceil(len(x) / float(scale)))
    rest = len(x) % scale
    div = len(x) / scale
    for i in range(2 * ns):
        if i == ns - 1:
            if rest != 0:
                print x[div * scale:]
            else:
                print x[(div - 1) * scale:]
        elif i == ns:
            if rest != 0:
                print x[:rest]
            else:
                print x[:scale]
        elif ns < i <= 2 * ns - 1:
            if rest != 0:
                print x[rest + scale * (i - ns - 1):rest + scale * (i - ns)]
            else:
                print x[rest + scale * (i - ns):rest + scale * (i - ns + 1)]
        elif 0 <= i < ns - 1:
            print x[i * scale:(i + 1) * scale]


if __name__ == '__main__':
    print "**************start**************\n\n"

    # expmt_div()
    # expmt_div(6, 3)
    # x = np.arange(0, 7.02, 0.01)
    # y = np.sin(x)
    # mfdfa = MFDFA(y)
    # mfdfa.setM(4)
    # mfdfa.process()
    # mfdfa.drawHq()
    # mfdfa.drawTq()
    # mfdfa.drawSpectrum()

    verify()
    # demo()

    print "\n\n**************end**************"
